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OSCILLATIONS OF MHD SHOCK WAVES ON THE SURFACES OF T TAURI STARS 
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ABSTRACT 

This work treats the matter deceleration in a magnetohydrodynamics radiative shock wave at the 
surface of a star. The problem is relevant to classical T Tauri stars where infalling matter is channeled 
along the star's magnetic field and stopped in the dense layers of photosphere. A significant new aspect 
of the present work is that the magnetic field has an arbitrary angle with respect to the normal to the 
star's surface. We consider the limit where the magnetic field at the surface of the star is not very strong 
in the sense that the inflow is super Alfvenic. In this limit the initial deceleration and heating of plasma 
(at the entrance to the cooling zone) occurs in a fast magnetohydrodynamic shock wave. To calculate 
the intensity of radiative losses we use "real" and "power-law" radiative functions. We determine the 
stability/instability of the radiative shock wave as a function of parameters of the incoming flow: velocity, 
strength of the magnetic field, and its inclination to the surface of the star. In a number of simulation 
runs with the "real" radiative function, we find a simple criterion for stability of the radiative shock 
wave. For a wide range of parameters, the periods of oscillation of the shock wave are of the order 
0.02 - 0.2 s. 
Subject headings: s 

tars: magnetic fields — stars: oscillations — MHD — accretion — shock waves — instabilities 



1. OVERVIEW OF THE PROBLEM 

In classical T Tauri stars (CTTSs) matter accretes from 
the disk to a star through magnetospheric funnel streams 
(Camenzind 1990; Konigl 1991; see also recent review 
by Bouvier et al. 2007). Similar type accretion but at 
smaller scales is expected to a magnetized white dwarf 
(e.g., Warner 1995) and magnetized neutron stars (e.g., 
Ghosh & Lamb 1979). In the funnel stream, matter is 
lifted above the equatorial plane and falls down onto the 
star due to the gravitational acceleration. Large-scale 
magnetospheric flow has been recently investigated in 2D 
magnetohydrodynamic (MHD) simulations (Romanova et 
al. 2002; Bessolaz et al. 2008) and in full 3D MHD simula- 
tions (Romanova et al. 2003, 2004; Kulkarni & Romanova 
2005). Many aspects of the global magnetospheric flow 
are now understood. However, interaction of the funnel 
streams with the surface of the star has not been ade- 
quately investigated. 

Theoretical models indicate that close to a star matter 
in the funnel streams is accelerated to almost free-fall ve- 
locity, before it hits the high-density layers of the stellar 
atmosphere. The matter rapidly slows down, forming a 
shock wave close to the stellar surface. Most of the en- 
ergy of the flow is radiated behind the shock wave (e.g., 
Lamzin 1995, 1998; Muzerolle et al. 1998; Calvet & Gull- 
bring 1998; Gullbring et al. 2000; Ardila & Basri 2000). 
In the case of CTTSs most of energy is radiated in the 
ultraviolet and soft X-ray bands (e.g., Calvet & Gullbring 



1998; Giinther & Schmitt 2008). The 3D MHD simula- 
tions of magnetospheric flow have show that the hot spots 
are inhomogeneous and are expected to have higher tem- 
perature in the central regions of spots compared to pe- 
ripheral regions (Romanova et al. 2004). This may have a 
number of important consequences for investigation of hot 
spots including the dependence of the filling factor on the 
wavelength. 

If star's magnetic field of the star is strong, then close 
to the stellar surface the magnetic energy-density is larger 
than that the kinetic energy-density of the matter. Con- 
sequently the matter is passively channeled along the field 
lines. In this sub-Alfvenic regime a hydrodynamic ap- 
proach is usually adopted for modeling the shock waves 
(see §2). The shock wave is found to be non-stationary. 
It oscillates with a high frequency due to the competition 
between accretion heating in the shock front and radiative 
cooling behind the front (Langer, Chanmugan & Shaviv 
1981; Chevalier & Imamura 1982). If the magnetic field 
is not very strong near the surface of the star then the 
flow may be super- Alfvenic. In this regime the orienta- 
tion of the magnetic field may influence stability of the 
shock. Only the special case of the flow perpendicular to 
the magnetic field has been considered so far. It is know 
that this transverse magnetic field can suppress instabil- 
ity of the shock wave (Smith 1989; Toth & Draine 1993). 
In this paper we investigate the stability of the radiative 
shock waves in the super- Alfvenic regime for different ori- 
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Fig. 1. — Sketch of the geometry of the flow. Matter with density p- ln and pressure p- ln flows towards the star with velocity 
v- m along the magnetic field which has strength B ln and is inclined relative to the normal to the surface of the star x at 
an angle \- Matter slows down in the MHD shock wave close to the surface of the star and radiates and cools down in 
the "cooling zone" . Matter may have different angle relative to the field, but we consider the coordinate system in which 
vectors velocity and magnetic field are parallel to each other. 



entations of the magnetic field relative to stellar surface. 
We consider small patch of the hot spot and investigate 
the stability of the radiative MHD shock wave for param- 
eters typical for CTTSs. Figure 1 shows a sketch of the 
considered geometry. 

For CTTSs the expected periods of the oscillations of 
the shock are very short. The periods vary between 0.02 
and 0.2 seconds depending on the parameters. Oscillations 
in this period range have not been observed so far. Smith, 
Jones & Clarke (1996) searched for rapid photometric vari- 
ability in several CTTSs in the range of periods from min- 
utes to hours and did not find oscillations. Much higher 
temporal resolution is required to resolve the oscillations 
discussed in this paper. 

Section 2 of this paper discusses the earlier research on 
radiative shocks. Section 3 discusses the model and ba- 
sic equations, and Section 4 the dimensionless variables 
and scalings. Section 5 comments on the dimensionless 
variables and the scalings of different quantities, and also 
describes the stationary structure of the shock. Section 6 
discusses the methods used to study the time-dependent 
shocks, and Section 7 gives our results. Section 8 gives the 
conclusions of this work. 

2. EARLIER RESEARCH OF RADIATIVE SHOCKS AND 
RADIATIVE COOLING FUNCTION 

The stability of shock waves has been investigated by 
different groups both analytically (linear analysis) and nu- 
merically. Most of the investigations have been restricted 
to a purely hydrodynamic analysis, because in many sit- 
uations matter is channeled along the field lines and the 
problem can be considered as non-magnetic. 

First results on interaction of the funnel streams with a 
star and cooling in the radiative shock wave have been ob- 



tained in application to accreting white dwarfs. Numerical 
modeling of radiative shock wave at the surface of white 
dwarf led to discovery of instability which is driven by 
alternation of accreting heating in the shock wave and ra- 
diative cooling behind the shock wave and in resulting os- 
cillations of the position of the shock front (Langer, Chan- 
mugan & Shaviv 1981). Linear analysis of the stability of a 
one-dimensional radiative shock wave was done by Cheva- 
lier & Imamura (1982). These authors assumed that the 
radiative cooling from a unit volume is p 2 T a , and they 
studied the stability as a function a. In a more general 
form the problem has been investigated by Ramachandran 

6 Smith (2004). This work assumed that the radiative 
losses vary as p^T a . The boundaries of the stable and 
unstable regions were found as well as the frequencies and 
growth rates of the lowest frequency modes for different 
values of a, (3. Ramachandran & Smith (2006) investi- 
gated the influence of the Mach number of the inflowing 
gas to the stability of the radiative shock wave. In ad- 
dition they considered flow at different adiabatic indices 

7 typical for astrophysical applications and two types of 
the boundary conditions at the "wall" where the flow is 
stopped. The influence of boundary conditions on the sta- 
bility of radiative shock waves has been investigated in 
detail by Saxton (2002). 

A detailed investigation of the linear and nonlinear evo- 
lution of one-dimensional radiative shock waves has been 
done by Mignone (2005, see also Toth & Draine 1993). 
The intensity of the radiative losses was taken to be p 2 T a . 
The stability has been investigated in the linear approxi- 
mation for the first eight low-frequency modes and it has 
been established, that: (1) The first eight modes are stable 
for a > 0.92; (2) The fundamental mode (n = 0) becomes 
unstable for a < 0.388, the n = 1 mode for a < 0.782; 
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(3) In the unstable regime the growth rate is larger for 
larger mode numbers, but the growth rate of the higher-n 
modes is not very different from the growth rate of the 
n = 7 mode; (4) The normalized frequencies of the corre- 
sponding modes have an approximately linear dependence 
on n, that is, oj n (a) = u>o(a) + nAw(a), and decrease as 
a increases (excluding the fundamental mode n = 0). 

The stability of a radiative shock wave in the presence of 
a magnetic field was investigated by Smith (1989). A more 
detailed investigation of the shock stability in presence of 
a magnetic field was done by Toth & Draine (1993; here- 
after TD93). For the situation considered, the matter flows 
onto the shock wave perpendicular to its front and the 
magnetic field has only component parallel to the front. It 
was concluded that even a modest magnetic field may lead 
to stabilization of a radiative shock wave which is unsta- 
ble in the hydrodynamic limit. The higher the harmonic 
number, the larger the value of the magnetic field which is 
needed for stabilization of the front for a fixed value of a. 
First of all, the magnetic field stabilizes the fundamental 
mode. In particular, for a — —0.5 the fundamental mode 
is stabilized at M^ 1 = [(B y / \J Air p) / v]i n = 0.15, while 
the higher modes are stabilized at — 0.5, where v 

is velocity of the incoming flow and B y is magnetic field 
parallel to the front. 

The subsequent investigation of the stability of the ra- 
diative shocks with transversal magnetic field has been 
done by Ramachandran & Smith, (2005; hereafter RS05). 
They considered different values of a, (3 in the dependence 
of the radiative losses, different values 7 for the adiabatic 
index, and also different Mach numbers in the inflowing 
matter. They obtained new results, and also rederived ac- 
curately results by TD93 for the case where 7 = 5/3 and 
= 2. 

These earlier studies show that the stability of radia- 
tive shock waves depends strongly on the functional de- 
pendence of the radiative cooling. In the present work 
we consider (3 = 2 and a monatomic gas with 7 = 5/3. 
We investigate the stability of the radiative shock waves 
using a "real" radiative loss function. We calculate this 
"real" function and approximate it with the power laws. 
We assume that there is collisional ionization equilibrium 
(CIE). In this approximation photons freely escape the 
plasma. Thus, full thermodynamic equilibrium is not es- 
tablished, and the Saha's formula (for calculation of the 
degree of the ionization) is not applicable (e.g., Spitzcr 
1968). The radiative losses under these conditions have 
been calculated by a number of authors. In these calcu- 
lations the abundances of the elements are assumed to be 
Solar. In the paper by Rosncr, Tucker, & Vaiana (1978; 
hereafter RTV) based on the calculations of Raymond & 
Smith (1977), the radiative losses in the temperature inter- 
val 10 4 3 K < T < 10 7 K are n e n H K(T), where n e ,n H are 
electron density and hydrogen density (total), while the ra- 
diation function A(T) has been approximated by a multi- 
segmented power-law (see below). Subsequently, A(T) was 
been calculated at the interval 10 3 65 K < T < 10 8 K (Peres 
et al, 1982). We refer to this function as the "real" cooling 



function and add a subscript RTV. The dependence is 
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The most recent results for radiative losses in the CIE- 
approximation are given by Gnat & Sternberg (2007; 
hereafter GS07). In this work it was accepted that 
the relative abundances of the hydrogen and helium 
are «Hc/ n H = 1/12. Abundances of other elements 
(C, Ni, O, N, Mg, Si, S, Fe) are small and they do not give a 
contribution to the total pressure. However, the intensity 
of the radiative losses significantly depends on the relative 
abundance of these elements. For the Solar abundance 
GS07 proposed an approximation formula: 



A GS = 2.3 x 10 



-19r 



-0.54 



erg cm s 



(1) 



in the temperature range 10 5 A' < T < 10 8 K. The left 
boundary of this interval corresponds approximately to 
the maximum of the "real" radiation function which is 
at 2.3 x 10 5 K. For temperatures higher than 6 x 10 7 K, the 
dominant mechanism of radiative losses is brcmsstrahlung 
radiation with the temperature dependence A ~ VT. 

According to calculations of GS07 in CIE- 
approximation, hydrogen is completely ionized for tem- 
perature T > 3 x 10 4 K, and Helium for T > 2 x 10 5 K. 
Accepting their abundances discussed above we obtain 
for T > 2 x 10 5 K, the electron density n e = «h + 2nne, 
and an average mass per particle is 0.6m p . Thus, the 
model which we use (based on the CIE-approximation) is 
applicable for T > 2 x 10 5 K. At lower temperatures, the 
partial ionization of Helium and the change of the average 
mass per particle become significant. At the temperature 
T = 3x 10 4 K the fractions of the neutral and ionized atoms 



of hydrogen and helium are: x(H) = 3.6 x 10 



-3 



s(H+) 



0.997, x(Re) = 0.4, a;(Hc+) = 0.6, a;(Hc++) = 0. Thus 
an average mass per particle is 0.62m p . At a temper- 
ature T = 2 x 10 4 K, the corresponding fractions are 
x(H) = 0.078, x(R+) = 0.922, z(He) = 0.993, a;(He+) = 
0.007, x(He ++ ) = 0, where an average mass per particle 
is 0.66m p . We neglect this factor. 

Subsequently we assume that plasma is an ideal gas with 
equation of state p = IZpT where 1Z = fcs/(0.6m p ) = 
1.385 x 10 8 erg g" 1 K _1 is gas constant. 

In the paper GS07 it is shown that for T < 10 6 K, the 
radiative losses, calculated in CIE-approximation, exceed 
losses which occur in non-stationary plasma cooling. This 
is connected with the fact that in the first case the ioniza- 
tion level of elements participating in the main radiative 
processes is lower because in the non-stationary regime re- 
combination lags the cooling. However, we use the more 
detailed approximation mentioned above after changing 
the normalization. 

Figure 2 shows the radiative cooling functions Ab.tv(T) 
(thick solid line) and Aqs(T) (thin solid line). The diag- 
onal dashed line shows the dependence of the upstream 
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Fig. 2. — Cooling function Artv(T') (thick solid line) and Ags(T") (thin solid line) as a function of T in K. The vertical 
dotted line shows our illustrative value of velocity before the shock v in cosx = 1-3 x 10 7 cm s^ 1 and the corresponding 
temperature behind the shock. 



velocity normal to the shock (v in cos x) on the temper- 
ature behind the shock (for 7 = 5/3, p m = 0): T s = 
3(w; n cosx) 2 /(167\L). The horizontal dashed line shows the 
velocity of the incoming flow v- m (for illustration of the 
calculated results we take «i n cosx = 1-3 x 10 7 cm s _1 ) 
and corresponding to this velocity temperature behind the 
shock wave front. 

Section 2 of the paper discusses the model and basic 
equations, and Section 3 the dimcnsionless variables and 
scalings. Section 3 comments on the dimensionless vari- 
ables and the scalings of different quantities. Section 4 
describes the stationary structure of the shock. Section 5 
discusses the methods used to study the time-dependent 
shocks, and Section 6 gives our results. Section 7 gives the 
conclusions of this work. 

3. FORMULATION OF THE PROBLEM 

We investigate formation and evolution of the shock 
wave near the surface of the star which forms as a result 
of disk accretion to a star through a funnel flow. The ac- 
creting matter is sufficiently ionized to satisfy the frozen-in 
condition so that matter of the funnel streams is channeled 
by the magnetic field and close to the star it flows along 
the field lines. 

Figure 1 shows the geometry. Matter with velocity Vi n , 
density p; n , and pressure p- ln flows towards the surface of 
the star along the magnetic field. The magnetic field has 
a strength B m and is directed at an angle \ relative to the 
normal vector to the surface of the star x. The x-axis is 
normal to the surface of the star and the y— axis is tangen- 
tial to its surface and is directed such that the magnetic 
field is located in the (x, y) plane. We neglect small per- 
turbations in z— direction associated with propagation of 
the Alfven waves. We consider that the magnetic field has 
an arbitrary inclination angle x relative to the normal to 
the star's surface. In general, the matter flow velocity is 
not parallel to the magnetic field. However, such paral- 
lel orientation can be obtained by transformation of the 



coordinate system. 

Heating of matter occurs in the front of the MHD shock 
wave. In the cooling zone behind the shock, matter radi- 
ates energy, is decelerated, and become denser up to the 
moment, when the radiative cooling stops. Formally this 
happens when T = 0. The height of the radiative zone is 
small compared to either width of the funnel stream or to 
the radius of the star. Thus we can neglect the inhomo- 
geneity of the accretion flow in the (y, z) directions. (See 
Canallc ct al. 2005 for a discussion of cases where the con- 
verging of the field lines is important.) We also neglect 
small effects associated with acceleration by gravity of the 
star (considered e.g., by Cropper et al. 1999) because the 
region we consider is small. Generally, the structure con- 
sisting of the shock wave and the cooling zone is unstable 
to both longitudinal, x) perturbations (which can be stud- 
ied in one-dimensional approach), and to transverse (y, z) 
perturbations (Bcrtschinger 1986; Imamura et al. 1996). 
The latter makes the problem more than one-dimensional 
even if the flow is homogeneous. In the present work the 
spatial perturbations arc not considered. 

We consider conditions where the star's mass is = 
0.8Mq and its radius is i?* = 2Rq. In this case the free- 
fall speed at the star's surface is Vff = y/2GM*/ 'R* = 
4 x 10 7 cm s _1 . If the temperature of accreting matter is 
10 4 K, then the sound speed is 2 x 10 6 cm s _1 (for aver- 
age mass per particle 0.6m p ). Clearly, the accretion is 
strongly supersonic; the sonic Mach number is 20. For 
a surface magnetic field B = 10 3 G and density of the 
inflowing matter 10 _11 g cm" 3 (Romanova et al. 2002), 
the Alfven velocity is ca = B/sJAirp = 10 8 cm s _1 . For 
these parameters the flow is sub-Alfvenic. However, for 
B < 3 x 10 2 G and the other parameters the same the 
flow is super-Alfvenic and perturbations from the shock 
wave cannot propagate up the stream. Strictly speaking 
we should compare the velocity of the accreting matter 
with the fast magnetosonic velocity, but it does not differ 
significantly from the Alfven velocity because the sound 
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Fig. 3.— The position of the boundary between stable and unstable radiative shock waves as a function of parameters 
(a, sinx) for a velocity of accretion V[ n cosx — 1-3 x 10 7 cm s _1 . 



speed is small. Thus, both cases are interesting: stability 
of the radiative shock waves for sub-Alfvenic and super- 
Alfvcnic inflows. However, we restrict the present study 
to super- Alfvenic inflow. 

We assume that the radiation zone behind the shock is 
optically thin in the direction perpendicular to the star's 
surface and to the front of the shock wave. The dominant 
radiation loss mechanisms are determined by the "ther- 
modynamic" state of plasma. We assume that the state 
of the matter can be described in terms of a temperature 
which is the same for the ions and electrons. For the con- 
sidered conditions, the radiative losses per unit volume is 
p 2 A(T), where p is the plasma density and A(T) is the ra- 
diative function obtained from Artv(T) or from Kqo,(T) 
by renormalization. We consider that the accreting matter 
is an ideal gas with adiabatic index 7 = 5/3 both before 
the compression at the shock front and in the radiation 
zone. 

In the absence of the magnetic field, and when the mat- 
ter falls perpendicular to the surface of the star, the tem- 
perature behind the shock is T — 10 b v 2 K, where v-j = 
w/(10 7 cm s" 1 ), and the average mass per particle is 0.6m p . 
For velocities of the incoming matter v = 4 x 10 7 cm s _1 , 
the temperature behind the shock reaches 1.6 x 10 6 K. At 
such temperature both hydrogen and helium atoms are 
completely ionized according to the CIE approximation. 

We consider the situation where the matter accretes to 
the star along magnetic field lines inclined by an angle \ 
to the normal to the star's surface. The shock is parallel 
to the star's surface and it remains parallel as it moves. 
That is, we consider variations only in the x— direction. 
The equations describing this situation are following: 



d(pv x ) 8 ( 2 Bl \ 



d{pVy) + I (pv x v y - ^ ] = 



dt dx 



Ait 



B x = const , 



013 d 

~~dt L + dx~ {VxBy ~ VyBx) = ' 



d ( pv 2 p B 2 



d_ 

dx 



dt V 2 7-1 8tt 



/ PV" IP \ By , 

Vx l~2~ + V^l ) + ^ X V ~ Vy ' 



(2) 

Here v 2 = v 2 x + v 2 . , B 2 = B 2 X + B 2 . In the unperturbed 



= VA(T) . 

, 2 = v 2 x + v 2 , B 2 = B 2 -<- R 2 
flow upstream of the shock, 



B x = B in cos x , By = B in sin x 



V x = -Win COS X , V y = -U in Sin \ , 



P = Pin , P= Ph 



dp d{pv x ) = 
dt dx 



We consider that the pressure in the incoming flow is very 
small so that it does not influence the MHD-shock wave 
and cooling zone. Equivalcntly, the sonic Mach number is 
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much larger than unity. As discussed earlier, we choose a 
coordinate system in which the flow velocity and magnetic 
field are both in the (x, y) plane. 

In contrast with the non-magnetic case (TD93) , the den- 
sity p does not increase unrestrictedly at the right-hand 
boundary of the radiative zone and correspondingly the 
velocity v x does not approach zero. Here, the tempera- 
ture T = p/(lZp) may approach zero not because p — > inf 
but because p — > 0. At the same time the total, gas plus 
magnetic field pressure remains constant. In the cooling 
zone the accreting matter does not reach zero flow speed. 
Thus, the radiative shock provides only part of the decel- 
eration and absorption of matter by the magnetized star. 

Our treatment of the stability of the radiative MHD 
shock wave is different in essential respects from that 
of TD93 and RS05. In these papers the magnetic field 
is transverse to the flow and parallel to the star's sur- 
face, B = B y y. We neglect the difference in cooling 
laws and adiabatic index 7 (RS05) which are not sig- 
nificant. To approach this case in our model, we need 
to take x ~> 7r /2, and parameters in the incoming flow 
v x = — v m cos x, By = Si n sin x — » £>in- The Alfvenic 
Mach number is 



M A 



Vin COS X 



cosx 



-B in /V47rp il 



(3) 



To have Ma a fixed value as x - * 7r /2, we let a = cosx/Ma 
and require 



B x = B in cos x -» , v x -> -M A 



B h 



y/4:irp in ' 



v y = —M A tanx 



Bu 



V^Kpi, 



00 . 



A tangential velocity v y can be included in the calculations 
of TD93 by a Galilean transformation to another reference 
frame. Thus, a radiative MHD shock wave with the mag- 
netic field parallel to the star's surface and perpendicular 
to the flow corresponds to the limit where x 7r /2 and 
a = cosx/Ma -> 0. 

4. DIMENSIONLESS VARIABLES AND SCALINGS 

Consider firstly the case where the radiation function is 
a power law, A(T) = A(lZT) a , where 1Z is the gas con- 
stant. The coefficient A has dimension cm 5 " a g _1 s Q_3 . 
For estimates we assume that the hydrogen and helium 
are completely ionized and that the average mass per par- 
ticle is 0.6m p . In the paper GS07 the radiative energy 
losses per unit of volume are given in the form of equation 
(1). With these simplifications and renormalization, 



Acs = 1.37 x lO^CfcT/cm 2 s~ z ) 



-°- 54 cm 5 



We rewrite our equations in dimensionless form choosing 
fiducial dimensional values for the main variables and in- 
troducing dimensionless variables A = A/Ao for different 
variables A. The fiducial values are taken to be: For the 
velocity, v = w; n cosx; for the density, p = p; n ; for the 
distance, v^ 2a /Apo; for time, v 2 T 2a / '(Apo); for the pres- 
sure, P0V 2 ; for the magnetic field, v ^/po; and for the tem- 
perature v 2 /TZ. In dimensionless variables, the system of 



equations has the same form. Subsequently we remove 
tilde signs from dimensionless variables. The formula for 
the radiative losses now has the form p 2 T a , while in the in- 
coming flow V in COS X — — 1, An = 1, B hl = (J\f^K, p in <C 1, 

where a is the inverse of the Alfven Mach number. 

It is clear that in power law case, A <~ T a , the 
stability/instability of the radiative shock in presence 
of a magnetic field does not depend solely on w in and 
B ln , but instead on their combination in the form 
(Bin/ Pin)/ Vin = c = M^ . Furthermore, the stability 
criterion does not depend on the coefficient A. However, 
A determines the spatial scale of the radiative zone and 
the temporal scale for oscillations if they are present. 

We chose the spatial scale is 10 7 v 4 - 08 /(/9_ii)cm and the 



time scale is 



/(/3_ii)s, where v-j is the velocity in 



units of 10 7 cm s _1 and p_n is the density in units of 
10~ n g cm~ 3 . In reality, both estimates are about three 
orders of magnitude larger than the observed values. Such 
scales follow from the solution of the equation for the sta- 
tionary shock wave, where the pressure goes to zero for 
x » 1 (TD93). The same is true for the scale of time. 
Thus, the height of the cooling zone is in fact ~ 10 4 cm, 
and the periods of oscillation are several hundredths of a 
second. 

We now consider the "real" radiative function (RTV) 
which is not a power law. This is why strictly speaking 
the stability condition of radiative shock waves depends 
not only on a but also on both v ln and i?i n . However, 
over a sufficiently wide range of temperatures, the radia- 
tive function can be roughly approximated by a power law 
of the temperature. This suggests that the stability con- 
dition will depend mainly on a with a weaker dependence 
on and B- m . 

5. STATIONARY STRUCTURE 

We are interested in the stability of the stationary flow 
consisting of the MHD shock wave and the downstream 
radiative zone. The stationary flow is described by 

PV X = ~j = -pin Win COSX , 



B B 2 

P + P v l + g^T = Qx = Pin + Pinvfn COS 2 X + ^ X , 



pV x V y 



B X By 

47T 



PinVi, 



B in 

4tt 



sin x cos x 



v x B. 



J 



dx 



VyB X = 



lp/p \ 

7-1 J 



p 2 A. 



(4) 



In the last equation there is no term describing the energy 
flux of the electromagnetic field, because the Poynting flux 
in our coordinate system is zero (E = v x B/c = 0). 

The stationary structure of the shock wave and radia- 
tive zone is described by the following equation for specific 
volume V , 



dV_ _ A(T) 
~dx ~ jV 2 F(V) 



(5) 
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N 


cooling 
function 


Vin COS X 

cm s -1 


sinx 


Bin 

(G) 


stab/ 
unstab 


I 


RTV 


1.3 x 10' 


0.154 


40.8 


unst 


II 


GS 


1.3 x 10' 


0.154 


40.8 


unst 


III 


RTV 


3 x 10' 


0.154 


94.9 


unst 


IV 


RTV 


1.3 x 10' 


0.368 


36.4 


stab 



Table 1 

The Table shows parameters used in our main runs. 



where 



F(V)=.fV 



1 + 



{B 2 jAn-j 2 Vf 



7+1 



q x - 2fV - 



qlB 2 x /8n 
{B 2 x /An-j 2 Vf 



+ 



(6) 



The other variables are determined from equations (4) . In 
particular, 



-JV , By 



B x q y 



Bl/An-fV 



By 

v v = -fT v x = 



B 2 /An-j 2 V ' 



P 



fv 



q 2 y B 2 /8n 



(B 2 /An-j 2 V) 2 



K v ; 



As mentioned we investigate stability of the structure 
"fast MHD shock wave" plus "cooling zone". Note that 
as the angle x between the flow velocity and the normal 
to the shock is decreased, the stationary MHD structure 
does not convert to the hydrodynamic one. In the limit 
X — > 0, q x i=a B 2 /(Ana 2 ). The pressure in the radiative 
zone is determined from the relation 



P 



Bj 
Ana 2 



-fv- 



q 2 y Bl, 



{B 2 j^-pvy 



(8) 



where q y oc sinx is small. For the considered conditions 
where a < 1, the denominator of the fraction of the last 
term goes to zero before the B 2 / '(Ana 2 ) — j 2 V goes to 
zero. Thus the pressure goes to zero for B 2 /An — j 2 V w 0. 
That is, the velocity of the flow at the exit from the ra- 
diative zone approaches the Alfven velocity and is v x w 
—B 2 /(Anj). Taking into account equation (7), we obtain 
from equation (8), 



Bl 



Bl Bt 



( v y 



Ana 2 An 8n \B 2 jAnj 



(9) 



^From this relation we obtain 



B„ 




— ^2(1-^) • (10) 
a 



Thus, as x — ► the components of the velocity and mag- 
netic field parallel to the shock approach finite values. In 
contrast, in the gas dynamic case these components ap- 
proach zero. The observed behavior is similar to that 
known to occur in MHD switch-on shock waves where fi- 
nite tangential velocity and magnetic field components are 
generated, and where the velocity of the flow behind the 
front is Alfvenic (Smith, 1993). Note that the parallel 
MHD shock wave becomes non-evolutionary, and it is re- 
placed by a switch-on shock wave for [v x /{B x /- s /Anp)\ ln < 
2 for 7 = 5/3 and p m = 0. 

6. METHOD 

We study the stability/instability of radiative shock 
waves in the presence of the magnetic field by integrat- 
ing the time-dependent one-dimensional MHD equations 
in a region containing the shock and the radiative cooling 
zone. We used an Eulerian variables. Consequently, the 
simulation region is chosen large enough to contain both 
the shock wave and the remote part of the radiative zone 
where energy losses are negligibly small. The location of 
the right-hand boundary of the simulation region is chosen 
so that the shock wave did not leave the region during os- 
cillations. For the calculations of the MHD flows we used a 
high resolution Godunov type numerical scheme (see, e.g., 
Kulikovskii, Pogorelov & Semenov 2000). 

For initial conditions we take the stationary flow with 
the shock wave and the radiation zone. At the top of 
the simulation region (see Figure 1) we determined either 
parameters of unperturbed flow (density, velocity, etc. ) 
if the incoming flow is super-Alfvenic, or we had "free" 
boundary conditions, if incoming flow is sub- Alfvenic. 

At the bottom of the simulation region (closer to stel- 
lar surface), we have a "layer" with no radiative losses. 
At this boundary we fixed the longitudinal velocity at a 
small value corresponding to the stationary solution. All 
other variables at the boundary had the same value as in 
the previous cell. Simulations have shown that results of 
modeling are not sensitive to the boundary conditions on 
the transverse components of velocity and magnetic field. 
At the boundary cell the density varied around some av- 
erage value and the mass did not accumulate there. For 
example in case I the sound speed in this cell has been 
(1 — 3) x 10 6 cm s _1 , the size of the grid = 81cm, the sound- 
crossing time < 10~ 4 s. We observed from simulations that 
variation of this boundary did not influence much to the 
oscillations. 

The spatial resolution has been chosen so that the radia- 
tion zone is covered by 200 cells. The size of the simulation 
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Fig. 4. — Time dependence of the shock front coordinate x s (top panel) and luminosity J (bottom panel) for case I 
("real" cooling function, v; n cosx = 1.3 x 10 7 cm s _1 ). The shock front coordinate is in units A, luminosity in units of 
J; n , and time is in seconds. 

region has been chosen such that during the oscillations the 
shock wave stayed inside the simulation region. Depend- 
ing on the amplitude of oscillations the simulation region 
incorporated from 500 to 1,500 cells. The time-step has 



been chosen automatically such that the Courant number 
is 0.5. 

7. RESULTS 



vo 
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Fig. 5. — Distribution of different parameters along the a;— axis: velocity component along the field lines (top panel), 
density and temperature (middle panel), magnetic field and plasma parameter /? (low panel). All variables are shown in 



dimensionlcss form. In the incoming flow v x 
from the right boundary, a star is at the left. 



-1, p = 1, T = (almost zero), B x = 0.48, B y = 0.076. Matter inflows 
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J(t) 




1.8.. 1.9 

time 



Fig. 6. — Time dependence of shock front coordinate x s (top panel) and luminosity J (bottom panel) for case II (power- 
law cooling function, Vj n cosx = 1.3 x 10 7 cm s _1 ). The coordinate of the shock front x s is in units A, the luminosity is 



in units of J; n , and the time is in seconds. 



We investigate the stability/instability of the radiative 
shock wave as a function of two main dimcnsionlcss param- 
eters: the inverse Alfven Mach number of the unperturbed 
upstream flow a = M7 1 and the inclination angle of the 
flow, Xi relative to the normal to the star's surface The 



limit (7 = corresponds to zero magnetic field. 

We performed a series of calculations in which we var- 
ied these two parameters. Figure 3 summarizes the re- 
sults. In the plane (a, sin \) markers show the parameters 
where calculations were done. The "squares" show the 




Fig. 7. — Time dependence of the shock front coordinate x s (top panel) and luminosity J (bottom panel) for case III 
("real" cooling function, v m cosx = 3 x 10 7 cm s _1 . The shock front coordinate is in units A, the luminosity is in units of 
Ji n , and the time is in seconds. 
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Fig. 8. — Fourier amplitudes of shock front velocity v s (top panel) and luminosity J (bottom panel) for case I. Velocity 
is in f; n cosx, luminosity is in J; n , frequency is in s _1 . 



parameters for which the radiative shock wave is stable 
and the "triangles" show cases where it is unstable. The 
solid straight line with a dashed continuation is the stabil- 
ity/instability boundary. The condition for stability can 
be expressed as 



3.7<7+ 1.4 sin x > 1 , 



(11) 



for conditions where the magnetic field is not very weak, 
a > 0.02. 

For weak magnetic fields (a < 0.02), the boundary 
of stability is located close to the vertical axis in Fig- 
ure 3. This part of the boundary is shown as a dashed 
line. Thus, in the absence of a magnetic field, the radia- 
tive shock wave (for the same parameters of the incoming 
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Fig. 9. — Fourier amplitudes of shock front velocity v s (top panel) and luminosity J for case II (cooling function AqS: 
Wincosx = 1-3 x 10 7 cms -1 )). The velocity is in units of fi n cosx, the luminosity is in units of J mi and frequency is in s _1 . 
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Fig. 10. — Fourier amplitudes of shock front velocity v s and luminosity J for case III (cooling function Artv 
fin cos x = 3 x 10 cms -1 )). The velocity is in units of v- m cos\, the luminosity is in units of J m , and the frequency 



is m s 



flow) is unstable. However, even a small magnetic field, 
in particular, inclined one, stabilizes the radiative shock 
wave. For a > 0.02 (Ma < 50), we find stability for an- 
gles x > 45°. If the magnetic field is sufficiently strong 
a > 0.27 (Ma < 3.7), then the radiative shock wave is 
stable for all angles \- 

We emphasize again that in the limit of small inclina- 
tion angles, \ 0: our model docs not correspond to a 
gas-dynamical flow with unperturbed magnetic field. 

To understand the temporal characteristics of unsta- 
ble radiative shock waves, we performed calculations us- 
ing the "real" cooling function Artv and for some cases 
with the function Aqs(T). In all calculations we followed 
the location of the shock wave x s (t) and the total inten- 
sity of radiation from the radiative zone per unit area 
J(t) = Jdxp 2 A. 

Figure 4 shows results of calculation of evolution of 
the radiative shock wave with radiation function Artv 



for the following parameters of the incoming flow: v- nl = 
1.3xl0 7 cm s"\ p hl = 10- n gcm- 3 , B in = 40.8G, sinx = 
0.154. Also, a = 0.13. The width of the cooling zone in the 
stationary regime is A w 1.64 x 10 4 cm, and the energy- 
density (excluding the much smaller thermal energy) is 
J in = (l/2)p in uf n cosx = 1-07 x 10 10 erg cnr 2 s -1 ). We 
point out once again that the incoming flow velocity is 
along the magnetic field, so that there is no Poynting flux. 

For the mentioned parameters, the stationary shock is 
unstable and the shock and radiative zone oscillate. The 
position of the front of the shock wave oscillates with a pe- 
riod w 0.025 s. The amplitude of oscillations of the front is 
modulated and varies with period s=s 0.3 s. The top panel 
of Figure 4 shows the position of the front as a function 
of time at the time-interval approximately equal to twice 
period of modulation. The bottom panel of Figure 4 shows 
the temporal variation of the radiation intensity from the 
radiation zone for the same time-interval as the top panel. 




0.4 .. 0.6 

time 



Fig. 11.— Time dependence of luminosity J for case IV (cooling function Artv , fin cos x 
luminosity is in units of Ji n , and the time is in seconds. 



1.3 x 10 7 cm s- 1 ). The 
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The position of the shock front is normalized to the width 
of the stationary cooling zone A. The intensity of radia- 
tion is normalized to the energy-density in the incoming 
flow Jj n . The time is in seconds. 

Figure 5 shows spatial distribution of the longitudi- 
nal (x— direction) velocity (top panel), density and tem- 
perature (middle panel), longitudinal magnetic field and 
plasma parameter (5 = 8irp/(Bl + By) (bottom panel) at a 
time corresponding to the maximum distance of the shock 
from the surface of the star. In stationary regime the shock 
wave is located at x = 0, while the radiative zone is below 
this. See animations at http://www.astro.cornell.edu/us- 
rus/shock.htm 

Figure 6 shows results of calculation of the evolution of 
the radiative shock wave for the same parameters in the 
incoming flow but with the radiation function Acs (case 

II) . The width of the cooling zone in the stationary regime 
is A w 2.65 x 10 4 cm. This figure shows position of the 
front of the shock wave x s (t) (top panel) and intensity of 
radiation J(t) (bottom panel). One can see that qualita- 
tively the oscillations are similar. However, the main fre- 
quency of oscillations is different: wrtv ~ 240s -1 , whereas 
lugs ~ 140s -1 . With the radiation function Acs, the mod- 
ulation of the oscillation amplitude of the position of the 
front is small. 

Figure 7 shows results of calculation of evolution of 
the radiative shock wave with radiation function Artv 
for the following parameters of the incoming flow (case 

III) : v in = 3 x 10 7 cm s" 1 , p in = 10 -11 gcm -3 , B in = 
94.9G, sinx = 0.154. Also, a = 0.13. The width of the 
cooling zone in the stationary regime is A w 8.86 x 10 5 cm, 
and the energy-density (excluding the small thermal en- 
ergy) is 1.35 x 10 11 erg cm -2 s -1 ). As in case I for these 
parameters the stationary shock is unstable and the shock 
oscillates. The position of the front of the shock wave 
oscillates with period « 0.21s. The top panel shows the 
position of the front as a function of time, and the bottom 
panel shows the variation of the radiation intensity for the 
same time interval as the top panel. 

In case III we increased both the velocity and magnetic 
field with the aim of checking our hypothesis, that the 
qualitative solution of the problem (while using the "real" 
radiation function Artv) is determined by the dimension- 
less parameters of the problem, sin \, and the magnetiza- 
tion, a, and not by dimensional values of velocity of the 
incoming flow and the magnetic field. 

Figure 8 shows results of the Fourier- analysis of the 
speed of the shock front and total luminosity per unit of 
area v s (u>) and J{uo). One can see that when modelling 
with radiation function Artv (case I) then two nearby 
maxima of the Fourier amplitudes are observed at the fre- 
quencies uj\ = 240s -1 , uj 2 — 260s -1 . We suggest that the 
combination of these two frequencies gives the amplitude 
modulation evident in Figure 4. 

Figure 9 shows the Fourier amplitudes for case II where 
we use Aqs- The Fourier-amplitudes for both, v s , and J 
have sharp maxima at frequencies divisible by the main 
frequency. 

Figure 10 shows the Fourier amplitudes for case III 
where we use Artv- The Fourier-amplitudes for both, 
v s , and J have sharp maxima at frequencies divisible by 
the main frequency The highest peak in Figures 8-10 cor- 



responds to the main oscillation frequency of the shock. 
The other peaks correspond to higher harmonics due the 
oscillations being anharmonic. 

We also investigated the stable regime of the radiative 
shock wave (case IV). To study the damping of the oscilla- 
tions we introduced a small (10%) perturbation of the ve- 
locity in the upstream flow during a limited time. We used 
the "real" radiative function, and parameters as shown in 
the Table for case IV. One can see that when perturbations 
reached the front of the shock wave, the shock wave be- 
gan to oscillate. However, these oscillations damped dur- 
ing several periods and stationary flow was re-established. 
Figure 11 shows an example of such damping. 

8. CONCLUSIONS 

This work has studied the stability/instability of the 
radiative MHD shock waves in the funnel streams of clas- 
sical T Tauri stars. Matter flowing to the surface of the 
star along the magnetic field is decelerated in the radiative 
shock wave. The shock may be stable or unstable. In the 
case of instability the shock position and other variables 
oscillate and this can give observable short time-scale vari- 
ability in the emitted radiation. A significant new aspect 
of the present work is that the magnetic field and the flow 
velocity parallel to it can have an arbitrary angle with 
respect to the normal to the star's surface. 

The shock wave has been modeled by solving the time- 
dependent MHD equations in one dimension (perpendicu- 
lar to the star's surface) taking into account the radiative 
losses. For the radiative losses we used either the "real" ra- 
diative function, approximated by segments of power laws 
(RTV) or by the power law function proposed by GS07. 

Results of modelling of the radiative shock waves show 
that there is a simple criterion of the shock stability: 
3.7a + 1.4 sinx > 1- This is for the case where the inflow 
to the shock is i>i„cosx = 1-3 x 10 7 cm s -1 . We believe 
that this criterion will not change significantly at larger 
inflow velocities. 

Comparison of the simulation results with the "real" 
(RTV) and power-law (GS) radiative functions shows that 
the qualitative results are similar. However, the periods of 
oscillations are significantly different. 

The periods of oscillations are of the order of hundredths 
of a second for w; n cosx = (1.3 — 3.0) x 10 7 cm s -1 . This 
period is expected to increase with V{ n . We estimate that 
P w 6 x 10 -3 [(wi n cosx)/(10 7 cm s -1 )] 3 s. The period of 
the oscillations varies, P = 0.02 — 0.2, depending on pa- 
rameters. 

Global three-dimensional simulations of magnetospheric 
accretion through the funnel streams have shown that hot 
spots on the surface of the star are not homogeneous: most 
of the kinetic energy flows in the central regions of the fun- 
nel stream so that the central regions of the spots are ex- 
pected to be hotter (and also denser) compared to periph- 
eral regions (Romanova et al. 2004; Kulkarni & Romanova 
2005). Romanova et al. (2004) have shown that the spots 
may have very small filling factor at highest density and 
temperature (less than 1%) and much larger filling factor 
at smaller densities and temperatures (see Fig. 3 of Ro- 
manova et al. 2004). This fact has been recently confirmed 
observationally by Giinther et al. (2007) who have shown 
that the filling factor in X-ray is smaller compared with 
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UV and optical bands which confirmed the inhomogeneity 
of the hot spots. Future research should be done for anal- 
ysis of the stability of the global shock wave which would 
cover a significant part of the hot spot (not a small part 
as usually considered including this paper). 

If the magnetic field near the star has a complex ge- 
ometry (e.g. Valenti & Johns-Krull 2004; Gregory et al. 
2006; Donati et al. 2007; Long, Romanova & Lovelace 
2007, 2008), then it is likely that some field lines arc in- 
clined to the surface of the star as considered in this pa- 
per. The present analysis is thus applicable to the sta- 
bility/instability of the shocks. If the complex field has 
significant transverse component then it may suppress os- 
cillations. 

Kravtsova & Lamzin (private communication) searched 
for oscillations of the shock in RW Aur using Crimean Ob- 
servatory facilities. They did not find oscillations, though 
their time-resolution was low (AP = 0.5 — Is. They con- 
cluded that in this star there are no oscillations with pe- 
riods P > 2s with amplitudes > 5% above the noise level. 
Higher time-resolution observations in a larger sample of 
CTTSs are needed to obtain a conclusive answer. It would 



be useful to have high (ms) time-resolution observations in 
the UV and X-ray bands in stars with high veiling, such as 
RW Aur and others, because these wavebands would cor- 
respond to oscillations of the central part of the hot spots 
(Romanova et al. 2004; Giinter et al. 2007) which may go 
into global oscillation mode with higher probability com- 
pared to peripheral parts observed in the optical band. A 
search in the optical band may bring interesting results as 
well because we do not know the details of the interaction 
of the funnel stream with a star on a global scale of the size 
of hot spot. To understand such physics, observations of 
variability time-scales would be informative and may help 
to shed a light to process of the funnel-star interaction. 
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